function [pars,Nu]                      = setup_function(setup_inputs,random_state,which_moments)

Nu                                      = [];

% Initialize random number generation
rng                                     default;
rng(random_state);


% Horizons for solultion, impulse response, and estimation, in months
pars.T                                  = 12*10;                                                                                % Date after which mean-reversion goes to zero
pars.tIRF                               = 8;                                                                                    % IRF horizon
pars.tSMM                               = 14;                                                                                   % Number of months for SMM analysis      

% Structural parameters - fixed in the analysis, except for one counterfactual
pars.r                                  = -1/12*log(0.70);                                                                      % Monthly model with annual discount rate of 0.70 
pars.M_c                                = 1;                                                                                    % Cash stock average            
if isfield(setup_inputs,'mean_reversion')
    pars.mean_reversion                     = setup_inputs.mean_reversion;                                                      % Shock persistence, expressed as the number of days to get back to within 10% of initial value
else
    pars.mean_reversion                     = 90;
end
pars.theta                              = -log(1-0.90)/( pars.mean_reversion/30 );                                              % Shock persistence, expessed as OU parameter                                            
pars.sigma_d                            = 0.010;                                                                                % Standard deviation of exposure to shock

% Structural parameters - estimated in the analysis
if isfield(setup_inputs,'C')
    pars.C                                  = setup_inputs.C;
else
    pars.C                                  = 0.063;
end
if isfield(setup_inputs,'S')
    pars.S                                  = setup_inputs.S;  
else
    pars.S                                  = 0.246;
end
if isfield(setup_inputs,'sigma')
    pars.sigma                              = setup_inputs.sigma;  
else
    pars.sigma                              = 0.039;
end
if isfield(setup_inputs,'M_e')
    pars.M_e                                = setup_inputs.M_e;  
else
    pars.M_e                                = 0.970;
end
if isfield(setup_inputs,'k')
    pars.k                                  = setup_inputs.k;
    pars.xi                                 = (1+pars.theta/(pars.r+pars.k))^(-1)*pars.sigma;
elseif (~isfield(setup_inputs,'k'))&&isfield(setup_inputs,'xi')
    pars.xi                                 = setup_inputs.xi;
    pars.k                                  = (pars.sigma/pars.xi-1)^(-1)*pars.theta - pars.r;                                     
else
    %pars.k                                  = 0.163;
    %pars.xi                                 = (1+pars.theta/(pars.r+pars.k))^(-1)*pars.sigma;
    pars.xi                                  = 0.00850;
    pars.k                                  = (pars.sigma/pars.xi-1)^(-1)*pars.theta - pars.r; 
end

% Check that parameters are consistent
if ~( abs( (1+pars.theta/(pars.r+pars.k))^(-1)*pars.sigma - pars.xi ) < eps )
    save('runs/setup_function_call.mat')
    error('Inconsistent parameters in setup_function')
end

% Parameters related to grids
Nu.NX                                   = 1e1+1;
if isfield(setup_inputs,'N')
    Nu.N                                    = setup_inputs.N;                                                                   % (1/2)*( # of gridpoints for M - 1 )
else
    Nu.N                                    = 15;
end
Nu.Deltat                               = 1/(Nu.N*pars.theta);                                                                  % Time period, expressed as fraction of a month
Nu.theta_target                         = pars.theta;

% For seeding random number generators
Nu.globalrandseed                       = 1;                                                                                    % For rng re-initialization at each call of produce_moments.m

% Parameters related to size of individual simulations
Nu.Nsim                                 = 512;                                                                                  % Number of initial states (# districts)
if isfield(setup_inputs,'shock_date')
    Nu.shock_date                           = setup_inputs.shock_date;                                                          % Day in the month at which the shock hits
elseif strcmp(which_moments,'baseline')
    Nu.shock_date                           = 8;
else
    Nu.shock_date                       = 25;
end
Nu.exposure                             = randn(Nu.Nsim,1);                                                                     % District-level exposure to aggregate shock (fixed across sims)
Nu.sampling_date                        = 30;                                                                                   % Sampling frequency, in days
if isfield(setup_inputs,'InitX')
    Nu.InitX                                = setup_inputs.InitX;
else
    if isfield(setup_inputs,'X0_mean')
        Nu.X0_mean                              = setup_inputs.X0_mean;
    elseif strcmp(which_moments,'baseline')
        Nu.X0_mean                              = 0.005;
    else
        Nu.X0_mean                              = 0.010;
    end
    Nu.XO_std                               = 0.050;                                                                            % Initial std of the distribution of X_0
    Nu.X0_a                                 = Nu.X0_mean*( Nu.X0_mean*(1-Nu.X0_mean)/(Nu.XO_std^2) -1 );                        % Shape parameter of beta distribution
    Nu.X0_b                                 = (1-Nu.X0_mean)/Nu.X0_mean*Nu.X0_a;                                                % Shape parameter of beta distribution 
    Nu.InitX                                = betarnd(Nu.X0_a,Nu.X0_b,Nu.Nsim,1);            
end
Nu.B                                    = 1e2;                                                                                  % Number of replications of the bootstrap
Nu.Nregions                             = Nu.Nsim;
Nu.Nperiods                             = pars.tIRF;
for b=1:Nu.B
    Nu.indices1(:,b)                    = randsample(Nu.Nregions,Nu.Nregions,true); 
    Nu.indices3(:,b)                    = randsample(Nu.Nperiods,Nu.Nperiods,true);
    Nu.indices4(:,b)                    = randsample(Nu.Nregions,Nu.Nregions,true);
end
rng(Nu.globalrandseed);                                                                                                         % Re-initialize random number generator

% Parameters related to SMM and solution methods
Nu.Nreps                                = 2e1;                                                                                  % Number of repetitions of the Monte-Carlo for SMM
Nu.Seeds                                = 1:1:Nu.Nreps;                                                                         % Seeding for replication
Nu.eps                                  = {};                                                                                   % Random shocks used for MC simulation with filtering
Nu.Nt                                   = ceil(pars.tSMM/Nu.Deltat);
for nr = 1:Nu.Nreps
        Nu.eps{nr}                              = randn(Nu.Nt+1,Nu.Nsim);
end
Nu.solutionmode                         = 0;                                                                                    % =0: compute only policy functions; ~0: compute everything
Nu.mc_method                            = 0;                                                                                    % MC sim method; =0: don't use transition matrix; =1: use dtmc; ~1: use inverse transform
Nu.tol_solve                            = 1e-5;                                                                                 % Tolerance for bounds
Nu.Nvec                                 = 9e0;                                                                                  % Grid size for search for initial points
Nu.Nstart                               = 5;                                                                                    % Number of initial points to try for each value of (sigma,xi)

% Parameters related to diagnostics
Nu.Nreps_diags                          = 2e2;                                                                                  % Number of repetitions of the Monte-Carlo for SMM
Nu.Seeds_diags                          = 1:1:Nu.Nreps_diags;                                                                   % Seeding for replication
Nu.eps_diags                            = {};                                                                                   % Random shocks used for MC simulation with filtering
for nr = 1:Nu.Nreps_diags
        Nu.eps_diags{nr}                        = randn(Nu.Nt+1,Nu.Nsim);
end
Nu.JacobianStepSize                     = 1e-4;                                                                                 % Step size in Jacobian for numerical differentiation                   
Nu.rh                                   = 0.01;                                                                                 % Relaxation parameter for bounds

% Other estimation options 
Nu.options_patternsearch                = optimoptions( 'patternsearch',...
                                                        'Display','off',...
                                                        'PlotFcn',[],...
                                                        'MeshTolerance',1e-6,...
                                                        'Cache','on',...
                                                        'MeshContractionFactor',0.1,...
                                                        'MeshExpansionFactor',10,...
                                                        'PollMethod','GPSPositiveBasisNp1',...
                                                        'TolFun',1e-6);     % Estimation options
Nu.options_fmincon                      = optimoptions( 'fmincon',...
                                                        'algorithm','active-set',...
                                                        'Display','off',...
                                                        'MaxFunEvals',2e4,...
                                                        'StepTolerance',1e-10,...
                                                        'TolFun',1e-10);     % Estimation options

% Moments
if nargin > 2 && strcmp(which_moments,'baseline')
    Nu.all_moments                          = {'SR','LR','IC','IC*LR','MSR','SRVAR','LRVAR','WNVAR'};
    Nu.targeted_moments                     = Nu.all_moments;
    Nu.Nmoments                             = numel(Nu.targeted_moments); 
    Nu.which_moments                        = 'baseline';
elseif nargin > 2 && strcmp(which_moments,'persist')
    Nu.all_moments                          = {'H1','H2','H3','H4','H5','H6','H7','H8'};
    Nu.targeted_moments                     = {'H2','H4','H8'};
    Nu.Nmoments                             = numel(Nu.targeted_moments); 
    Nu.which_moments                        = 'persist';
else
    Nu.all_moments                          = {'Hm5','Hm4','Hm3','Hm2','Hm1','H0','H1','H2','H3','H4','H5','H6','H7','H8'};
    %Nu.targeted_moments                     = {'Hm3','Hm2','Hm1','H0','H2','H4','H8'};
    Nu.targeted_moments                     = Nu.all_moments;
    Nu.Nmoments                             = numel(Nu.targeted_moments); 
    Nu.which_moments                        = '';
end
for m=1:Nu.Nmoments
        Nu.idx_targeted_moments(m)                                      = find(strcmp(Nu.all_moments, Nu.targeted_moments{m}));
end 
Nu.log_targeted_moments                 = false(numel(Nu.all_moments),1); 
for m=1:Nu.Nmoments                                                                
        Nu.log_targeted_moments(Nu.idx_targeted_moments(m))             = true;
end   

% Which parameters to estimate
Nu.all_pars                             = fieldnames(pars);
Nu.full_pars                            = {'S','xi','sigma','M_e','C'};
if isfield(setup_inputs,'targeted_pars')
    Nu.targeted_pars                        = setup_inputs.targeted_pars;
else
    Nu.targeted_pars                        = {'S','M_e','C'};
end
Nu.Npars                                = numel(Nu.targeted_pars);
Nu.fixed_pars                           = setdiff(Nu.all_pars,Nu.targeted_pars);

% Joint constraints on all 5 parameters
Nu.lb_k                                 = 0.05;
Nu.ub_k                                 = min( 1/(Nu.Deltat*Nu.NX), (1/2)*(1/Nu.Deltat - pars.r) );
if isfield(setup_inputs,'lb_sigma')
    Nu.lb_sigma                             = setup_inputs.lb_sigma;
elseif strcmp(which_moments,'baseline')
    Nu.lb_sigma                             = 0.020;
else
    Nu.lb_sigma                             = 0.010;
end
if isfield(setup_inputs,'ub_sigma')
    Nu.ub_sigma                             = setup_inputs.ub_sigma;
elseif strcmp(which_moments,'baseline')
    Nu.ub_sigma                             = 0.080;
else
    Nu.ub_sigma                             = 0.100;
end
Nu.lb_xi                                = max(0.0050,(pars.r+Nu.lb_k)/(pars.r+Nu.lb_k+pars.theta)*Nu.lb_sigma);
Nu.ub_xi                                = min(0.0400,(pars.r+Nu.ub_k)/(pars.r+Nu.ub_k+pars.theta)*Nu.ub_sigma);
if strcmp(which_moments,'baseline')
    Nu.lb_global                            = [0.15 Nu.lb_xi Nu.lb_sigma 0.95 0.04]';
else
    Nu.lb_global                            = [0.01 Nu.lb_xi Nu.lb_sigma 0.92 0.04]';
end
if isfield(setup_inputs,'C')
    if setup_inputs.C < Nu.lb_global(end)
    Nu.lb_global(end)                       = 0;
    end
end
if isfield(setup_inputs,'S')
    if setup_inputs.S < Nu.lb_global(1)
    Nu.lb_global(1)                         = setup_inputs.S;
    end
end
if strcmp(which_moments,'baseline')
    Nu.ub_global                            = [0.30 Nu.ub_xi Nu.ub_sigma 0.99 0.08]';
else
    Nu.ub_global                            = [0.40 Nu.ub_xi Nu.ub_sigma 0.99 0.08]';
end
Nu.A_global                             = [0         -sqrt(Nu.N/pars.theta)             0                                               1              1;
                                           0         -sqrt(Nu.N/pars.theta)             0                                              -1              0;
                                           0         -1                                 (pars.r+Nu.lb_k)/(pars.r+Nu.lb_k+pars.theta)    0              0;
                                           0          1                                -(pars.r+Nu.ub_k)/(pars.r+Nu.ub_k+pars.theta)    0              0;
                                           0         -sqrt(Nu.N/pars.theta)             0                                               0              1/2];
Nu.b_global                             = [ pars.M_c;
                                           -pars.M_c;
                                           0;
                                           0;
                                           0];  

Xi                                      = [pars.S; pars.xi; pars.sigma; pars.M_e; pars.C];
if any(Nu.A_global*Xi > Nu.b_global)||any( Xi > Nu.ub_global )||any( Xi < Nu.lb_global )
        save('../runs/setup_call.mat');
        error('Error: point violates the global constraints \n');
end

% Constraints on the 3 estimated parameters
if isequal(Nu.targeted_pars,{'S','M_e','C'})
Nu.lb                                   = [Nu.lb_global(1)      max(Nu.lb_global(4),pars.M_c - sqrt(Nu.N/pars.theta)*pars.xi)   Nu.lb_global(5)]';
Nu.ub                                   = [Nu.ub_global(1)      Nu.ub_global(4)                                                 min(2*sqrt(Nu.N/pars.theta)*pars.xi,Nu.ub_global(5))]';
Nu.A                                    = [0    1    1];
Nu.b                                    = pars.M_c + sqrt(Nu.N/pars.theta)*pars.xi;
Xi                                      = [pars.S; pars.M_e; pars.C];
    if any(Nu.A*Xi > Nu.b)||any( Xi > Nu.ub )||any( Xi < Nu.lb )
            save('../runs/setup_call.mat');
            error('Error: point violates the local constraints \n');
    end
else
    if isequal(Nu.targeted_pars,{'S','sigma','M_e','C'})
        Nu.lb                                   = [Nu.lb_global(1)      Nu.lb_global(3) max(Nu.lb_global(4),pars.M_c - sqrt(Nu.N/pars.theta)*pars.xi)   Nu.lb_global(5)]';
        Nu.ub                                   = [Nu.ub_global(1)      Nu.ub_global(3) Nu.ub_global(4)                                                 min(2*sqrt(Nu.N/pars.theta)*pars.xi,Nu.ub_global(5))]';
        Nu.A                                    = [0    0    1    1];
        Nu.b                                    = pars.M_c + sqrt(Nu.N/pars.theta)*pars.xi;
        Xi                                      = [pars.S; pars.sigma; pars.M_e; pars.C];
        if any(Nu.A*Xi > Nu.b)||any( Xi > Nu.ub )||any( Xi < Nu.lb )
                save('../runs/setup_call.mat');
                error('Error: point violates the local constraints \n');
        end
    else
        if isequal(Nu.targeted_pars,{'S','sigma','M_e','C'})
            Nu.lb                                   = [Nu.lb_global(1)      Nu.lb_global(3) max(Nu.lb_global(4),pars.M_c - sqrt(Nu.N/pars.theta)*pars.xi)   Nu.lb_global(5)]';
            Nu.ub                                   = [Nu.ub_global(1)      Nu.ub_global(3) Nu.ub_global(4)                                                 min(2*sqrt(Nu.N/pars.theta)*pars.xi,Nu.ub_global(5))]';
            Nu.A                                    = [0    0    1    1];
            Nu.b                                    = pars.M_c + sqrt(Nu.N/pars.theta)*pars.xi;
            Xi                                      = [pars.S; pars.sigma; pars.M_e; pars.C];
            if any(Nu.A*Xi > Nu.b)||any( Xi > Nu.ub )||any( Xi < Nu.lb )
                    save('../runs/setup_call.mat');
                    error('Error: point violates the local constraints \n');
            end
        else
            if isequal(Nu.targeted_pars,{'S','xi','sigma','M_e','C'})
                Nu.lb                                   = Nu.lb_global;
                Nu.ub                                   = Nu.ub_global;
                Nu.A                                    = Nu.A_global;
                Nu.b                                    = Nu.b_global;
            end

        end
    end     
end

end
